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ABSTRACT 

We calculate the cosmological evolution of the 1-point probability distribution function (PDF), 
using an analytic approximation that combines gravitational perturbation theory with the Edge- 
worth expansion of the PDF. Our method applies directly to a smoothed mass density field or to 
the divergence of a smoothed peculiar velocity field, provided that r.m.s. fluctuations are small com- 
pared to unity on the smoothing scale, and that the primordial fluctuations that seed the growth of 
structure are Gaussian. We use this "Edgeworth approximation" to compute the evolution of {6\6\) 
and these measures are similar to the skewness and kurtosis of the density field, but they are 
less sensitive to the tails of the probability distribution, so they may be more accurately estimated 
from surveys of limited volume. We compare our analytic calculations to the results of cosmological 
A^-body simulations in order to assess their range of validity. The Edgeworth approximation for 
the PDF, and the computations of {S\6\) and {\5\) that are based on it, remain quite accurate until 
the r.m.s. density fluctuation is ci ~ 1/2, or, more generally, until the magnitude of the skewness 
approaches one. The skewness and kurtosis of the density field stay remarkably close to the values 
predicted by perturbation theory even when the r.m.s. fluctuation is cr = 2. When cr <C 1, the 
numerical simulations and perturbation theory agree precisely, demonstrating that the A^-body 
method can yield accurate results in the regime of weakly non-linear clustering. We show analyti- 
cally that "biased" galaxy formation preserves the relation (5^) oc (5^)^ predicted by second-order 
perturbation theory, provided that the galaxy density is a local function of the underlying mass 
density. The constant of proportionality depends on the shape of the biasing function that relates 
galaxy and mass densities. Our results should be useful in the analysis of large-scale density and 
velocity flelds, allowing one to derive constraints on the nature of primordial fluctuations, the value 
of the cosmological density parameter, and the physical processes that govern galaxy formation. 

Subject Headings: cosmology: theory — large-scale structure of the Universe 
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1. Introduction 

Most theories for the formation of structure in the universe assume that this structure de- 
veloped by gravitational instability from small-amplitude primordial fluctuations. The simplest 
hypothesis is that these fluctuations were Gaussian, and simple versions of inflationary cosmology 
naturally produce Gaussian fluctuations from the quantum fluctuations of the inflaton field (Guth 
& Pi 1982; Hawking 1982; Starobinsky 1982; Bardeen, Steinhardt & Turner 1983). Topological de- 
fect models predict non-Gaussian fluctuations, as do some specialized versions of inflation involving 
multiple scalar fields (e.g. Zel'dovich 1980; Turok 1989; Barriola &: Vilenkin 1989; Allen, Grinstein 
& Wise 1987; Kofman &; Pogosyan 1988; Salopek, Bond Sz Bardeen 1989). Observational evidence 
for or against Gaussian primordial fluctuations can therefore provide important clues to physics in 
the very early universe, and to the physical origin of today's large-scale structure. 

In the linear approximation for gravitational instability, fluctuations that are initially Gaussian 
remain Gaussian. However, non-linear effects quickly distort Gaussian fluctuations, and they are 
quite significant on all scales that can be probed by current observational surveys. In this paper we 
show how to compute the evolution of the 1-point probability distribution function (PDF) in the 
weakly non-linear regime, incorporating systematically the first departures from Gaussianity. We 
also compute the evolution of two quantities that measure the asymmetry and width of the PDF, 
respectively {6\5\) and {\S\), where S = {p — {p))/{p) is the density contrast and the brackets (. . .) 
denote averaging over the PDF. These measures are similar to the skewness and kurtosis of the 
density field, but they may be less subject to observational sampling errors because they are less 
dependent on the high-(5 tail of the PDF. We present similar calculations for the divergence of the 
velocity field, and we consider the effects of "biased" galaxy formation on the moments and PDF 
of the density field. Throughout the paper wc compare our analytic calculations to cosmological 
A'^-body simulations in order to assess their range of validity. Our results should be useful in the 
analysis of large-scale density and velocity fields, providing tools with which to test the hypothesis 
of Gaussian initial fluctuations and constrain the value of $7, the cosmological density parameter. 

There have been previous efforts to compute the evolution of the PDF from Gaussian initial 
conditions, employing the Zel'dovich approximation (Kofman et al. 1993) as well as rigorous pertur- 
bation theory (Bernardcau 1992). However, these methods compute the probability distribution of 
the unsmoothed flnal density field that evolves from smoothed initial conditions. For comparison to 
observational data, the relevant quantity is the PDF of the smoothed final density field that evolves 
from unsmoothed initial conditions, which may be far into the non-linear regime on scales below the 
smoothing length. Because dynamical evolution is non-linear, the effects of gravity and smoothing 
do not commute. Padmanabhan & Subramanian (1993), recognizing this problem, have attempted 
to compute the PDF of the smoothed final density field via the Zel'dovich approximation, i.e. by 
using first-order perturbation theory in Lagrangian space. In a series of recent papers, we have 
shown how to calculate moments of the PDF of a smoothed final density field using perturbation 
theory; a rigorous perturbative calculation of order-n moments requires order- (n — 1) perturbation 
theory (Juszkiewicz & Bouchet 1992; Bouchet et al. 1992a; Juszkiewicz, Bouchet &: Colombi 1993; 
see also Goroff et al. 1986). In this paper wc combine these results with the Edgeworth expansion 
(see Cramer 1946) to obtain an approximation to the full PDF. 

2. The Edgeworth Expansion 

We wish to examine how gravitational instability drives the PDF away from its initial state, 
which we assume to be Gaussian. We first introduce the Gram-Charlier expansion, which allows 
one to reconstruct the PDF from its moments. We then summarize the predictions of perturbation 
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theory for the moments of the smoothed mass density contrast 5. Finally, we rearrange the Gram- 
Charlier series by collecting all terms of the same order. The result is the proper asymptotic 
expansion of the PDF in powers of cr, the standard deviation of 6. 

Our object of study is p{v), the PDF of the density field in terms of the standardized random 
variable v = bja. The probability that the density contrast at a randomly chosen location lies in the 
range v < bja < v -\- dv \s, 'p(v)dv. Let us also introduce <^{y) = {2tt)~^^^ exp(— 1/^/2), a Gaussian 
(or Normal) PDF. Since we want to describe evolution from Gaussian initial conditions, it makes 
sense to consider an expansion of in terms of and its derivatives. The Gram-Charlier 
series (Cramer 1946 and references therein) provides such an expansion: 

P{U) = CO m + ^ (Z.) + I 0(2) (j.) + . . . , (1) 

where q are constant coefficients. Superscripts denote derivatives with respect to v: 



<l^^'H^)^^ = {-^YHe{^)H^), (2) 

where Hi is the Hermite polynomial of degree i. Table 1 gives expressions for the first seven H£. 
The Hermite polynomials satisfy orthogonality relations (e.g. Abramowitz & Stegun 1964): 

£h,M;^„M = (3) 

Therefore, multiplying both sides of equation (1) by and integrating term by term yields 

/oo 
Hi{u)p{u)dv. (4) 
-oo 

Equation (4) gives cq = 1, ci = C2 = 0, while for the next four coefficients in the series we obtain 
Q = (-l)^5^c7^-2, for3<^<5; Cg = Sga^ + Sfa^ . (5) 

Here 

S^ = Me/a^'-\ (6) 
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with Mi being the ^:th cumulant^ of p{S). Each cumulant is a combination of the central moments 

/oo 
p{iy) dv . (7) 
-oo 



In particular, 



M2 = /X2 = g\ M3 = //3, M4 = - 3C7^ 

M5 = /X5 - lOa^Ma, M6 = H6- 15a^M4 - 10M| - 15(7^ . 



For a zero-mean, Gaussian distribution, all cumulants vanish except M2. Ratios of the cumulants to 
the appropriate powers of a provide convenient measures of deviations from the Gaussian shape. For 

example, the skewness M^/a^ measures the asymmetry of the distribution, while the kurtosis M^/a'^ 
measures the flattening of the tails relative to a Gaussian. The significance of the "normalized 
cumulants" will become clear shortly. 

Thus far our discussion has been entirely general. Now let us consider the application to the 
gravitational evolution of initially Gaussian fluctuations. In the weakly non-linear regime, when 
the r.m.s. amplitude of density fluctuations, a, is smaller than unity, the time evolution equations 
for the smoothed density contrast can be solved approximately by using the perturbative expansion 

6 = 61 + 62 + .... (9) 

Here 61 = 0{a) is the solution of the equations with all non-linear terms set to zero; 62 = 0{a^) is 
the solution of equations of motion with quadratic terms included iteratively by using 61 as a source, 
and so on (see, e.g., §18 in Peebles 1980). For an Einstein-de Sitter = 1) cosmology, all terms in 
this expansion are known (Fry 1984); for Q ^ 1, solutions are available only for the first few terms 
(Bouchet et al. 1992a, hereafter BJCP). We assume that 61 is a Gaussian field. All of its statistical 
properties are therefore determined by its power spectrum, P{k) = (5^), where k is the comoving 
wavcnumbcr. Since one constructs the higher order terms iteratively out of 5i, the initial power 
spectrum also determines the non-linear, dynamically driven deviations from Gaussian behavior. 
One can compute the moments of the density field needed for the Gram-Charlier expansion of p{v) 
by raising both sides of equation (9) to the appropriate power, then averaging. The cumulants 
Ml can then be expressed in terms of the parameters, multiplied by the appropriate powers 
of £7 = {6\)^/'^. For details of the calculations, see Goroff et al. (1986, hereafter GGRW) and 
Juszkiewicz et al. (1993, hereafter JBC). The lowest non- vanishing contribution to the skewness is 

Ms = ?,{5l62) = 0{cT^) , (10) 

and the £:th cumulant is of order (Fry 1984; Bernardeau 1992) 

M, = 0{a^'-^) . (11) 

Equation (11) implies that, in perturbation theory, the normalized cumulants S(, defined by equation 
(6) are always of order unity. In an Einstein-de Sitter universe, the Si remain independent of time 
in the perturbative regime, determined only by the slope of the power spectrum near the smoothing 
scale (GGRW; JBC). For J7 7^ 1 models, only ^3 has been calculated. Although ^3 does change with 
time in this case, the time dependence is extremely weak, and for identical spectra and smoothing 



^ A cumulant of arbitrary order I is given by d^ln(e*'') /dt^, evaluated at t = 0. Cumulants are 
often called "reduced" or "connected" moments. 
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Figure 1 — The derivatives of the Gaussian function that appear in the low-order terms of the Edgeworth series, 
equation (12). Sohd, dotted, and dashed lines show 4>'^^^(i'), and ^(^'(f) respectively, multiplied by factors 

that appear in equation (12). 

functions S3 remains within several per cent of the corresponding Einstein-de Sitter value provided 
0.05 <n<3 (BJCP; JBC). 

The normalized cumulants Se have both a dynamic and a static application: they describe 
the time evolution of moments of the PDF at a fixed smoothing scale, but they also describe the 
relation between moments of the PDF at a fixed time on different smoothing scales. In the latter 
case, one must also include the scale-dependence of the Se if the initial power spectrum is not 
scale-free. 

For our immediate purpose here, the important consequence of the fact that the normalized 

cumulants arc 0(1) in perturbation theory is that the Gram-Charlier series is not a proper asymp- 
totic expansion for p{i'). In an asymptotic expansion, the remainder term should be of higher order 
than the last term retained. However, if we truncated the series (1) at the (p^"^^ term, which is 
0(a"^), we would miss another 0(<t^) contribution coming from cq (equation [5]). In order to deal 
with this problem, let us rearrange the Gram-Charlier expansion by collecting all terms with the 
same powers of a. The result is the so-called Edgeworth series, with the first few terms given by 

pH = - ^53</.(^)(^^)^ + ^S,cf>^'^\u)a^ + ^5|</,(6)(^.)a2 + 0{a') . (12) 

Figure 1 plots the derivatives 4>^^\ cf>^^\ and (p^^^ that appear in equation (12). Cramer (1946) lists 
the Edgeworth series to higher order, and he proves that it is a proper asymptotic expansion. This 
proof is directly relevant to our purposes, since it implies that there are no additional O(cr^) terms 
hiding in the Gram-Charlier series at ^ > 6. 

Now we can see the attractiveness of the Edgeworth series for describing the gravitational 
evolution of Gaussian fluctuations: it becomes a scries expansion for the evolving PDF in powers 
of the r.m.s. fluctuation a. This makes physical sense because the Edgeworth series provides 
an expansion about a Gaussian probability distribution. If the initial fluctuations are Gaussian, 
then we expect the terms describing successively larger departures from a Gaussian PDF to come 
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in with successively higher powers of a. For similar reasons the Gauss-Hermite series, which is 
closely related to the Edgeworth expansion, has recently found applications in stellar dynamics as 
a description of galaxy line profiles (e.g. van dcr Marel & Franx 1993; Gerhard 1993). 

Given equation (12), we can compute the Edgeworth approximation to the PDF provided that 
we can compute Sf, to the required order. In this paper we will make use of the second-order 
approximation, 

1 + l-S^aH^{v) 



p{v) 

and the third-order approximation. 
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<j>{v), (13) 



4>{v). (14) 



Although equation (13) contains only a single explicit power of a, it is appropriately described 
as a second-order approximation because the parameter 5*3 remains zero until second order in 
perturbation theory. Similarly, equation (14) is a third-order approximation because the calculation 
of ^4 requires third-order theory. GGRW compute values of S'3, and ^5 for a cold dark matter 
(CDM) power spectrum smoothed with a Gaussian filter, and JBC derive 6*3 for power law spectra 
smoothed with Gaussian or top hat filters. Bernardeau (in preparation) derives for a power law 
spectrum and top hat filter. Note that all of these values are for the smoothed final density fields, 
so by combining them with the Edgeworth expansion one incorporates the effects of gravitational 
evolution and smoothing in the correct order. 

Thus far we have considered the PDF of the mass density field, which can be estimated from 
a galaxy redshift survey if one assumes that galaxies trace mass, or if one incorporates an assumed 
model of "biased" galaxy formation to describe the relation between the galaxy distribution and 
the mass distribution (see further discussion in §5). However, if the input data set is a peculiar 
velocity field, it makes more sense to look at the velocity divergence, 

e = V-v/Ho, (15) 

where v is the peculiar velocity and Hq is the Hubble constant. Our results above can be im- 
mediately generalized; one just adopts the values of Se that are relevant for 6 instead of S. The 
computation of these moment ratios is discussed by Bernardeau et al. (in preparation, hereafter 
BJDB), who give values of S'361 for power law spectra smoothed with Gaussian or top hat filters, 
and by Bernardeau (in preparation), who gives S^o for power law spectra smoothed with a top hat 
filter. 



3. Measures of Asymmetry and Width 

In Gaussian models, second-order perturbation theory predicts that = (<^'^)/(<^^)^ should be 
a constant, depending only on the shape of the power spectrum near the smoothing scale. One can 
check whether the density fields derived from galaxy redshift surveys and peculiar velocity surveys 
obey this relation in order to test the hypothesis of Gaussian initial fluctuations (for studies along 
these lines, see Bouchet et al. 1992b, 1993; Park 1991; Silk & Juszkiewicz 1991; Coles & Frenk 
1992). A disadvantage of using the third moment {5'^) is that it is quite sensitive to the tails of the 
PDF, so it is subject to sampling errors unless the survey volume is very large. Recognizing this 
problem, Nusser & Dekel (1993) use the quantity {5\8\) instead of {5^) as a measure of asymmetry in 
the PDF. In the Appendix, we outline a direct perturbative calculation of {6\6\) for the unsmoothed 
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density field that evolves from smoothed initial conditions. However, once one has computed the 
value of 5*3 from perturbation theory, one can derive the same result much more simply using the 
Edgeworth approximation to the PDF. Indeed, to calculate the first non-vanishing term in the 
asymptotic expansion for {S\S\), it is sufficient to use the second-order Edgeworth expansion (13): 

{d\d\) = \p{u) - p{-u)] diy=^ i^^H^iu) <t>{u) du + 0{a') , (16) 

with H2,{h') = — Su. Equation (16) is accurate to 0{a^) because the additional terms in the 
next order of the Edgeworth expansion are symmetric in u and make no contribution to 
Evaluating the simple integral on the right hand side of equation (16), we find 

m) = {^S,a^ + O(a^) . (17) 

This result is more general than the value for derived in the Appendix "from first principles" 

because it is not restricted to the unsmoothed field. The effects of a low pass filter can be included 
simply by using the value of S^, calculated from perturbation theory for the smoothed final density 
contrast. One can also incorporate the effects of shot noise, since the impact of Poisson fluctuations 
on ^3 can be calculated easily. Redshift-space distortions can be treated using the S-i results of 
Bouchet et al. (in preparation), and biased galaxy formation can be included using the results 
in §5 below (see also Fry & Gaztahaga 1993). Finally, from the above derivation it is clear how 
to compute the moment {0\9\) of the smoothed velocity divergence; one simply substitutes the 
appropriate value of S^e for ^3. 

The quantity {5\5\) measures asymmetry of the PDF, in similar fashion to the third moment 
[8'^). A measure analogous to the reduced fourth moment, M4 = {S"^) — Su^, is the expectation 
value of minus the contribution expected for a Gaussian PDF. Nusser & Dekel (1993) use this 
quantity, {\8\) - {2/tiY/'^(j, to measure the width of the PDF relative to a Gaussian distribution. 
We can compute its evolution in perturbation theory using the third-order Edgeworth expansion 
(14): 



/"OO 

{\5\)=CT / u\p{v)+p{-y)]du 

Jo 



2a 



/ i>(l){i>)dv + —S^a^ I i'Hi{i')(l>{i>)dv + —S^a^ I vHQ{v)(l>{i')dv 
Jo 4! Jo 6! Jo 



(18) 



Evaluating the integrals (the last two by parts, using equation [2]) and rearranging terms we find 

m - (2/7r)VV = {2n)-'/'{Si - S,)^ + 0{a') . (19) 

Once again, we can compute this quantity for the smoothed final density field or the smoothed 
velocity divergence by inserting the appropriate values of S3 and 84^ from perturbation theory. The 
fourth moment weights the tails of the PDF heavily, but {\5\) responds primarily to the width of 
the PDF near its peak, so a distribution with high kurtosis (high ^4) tends to have a negative value 
of {\5\) - (2/7r)VV. 

For convenience in the sections that follow, we now introduce the notation 




M3 = {m), Ss = ^ — Ss, (20) 

M4 = m) - (2/7r)^/V, 54 = (27r)-^/^(5| - S4)/12, (21) 
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in obvious analogy to the cumulants M3, M4 and the normalized cumulants S3, S4. In this notation, 
the predictions of perturbation theory combined with the Edgeworth approximation are simply 



M3 = Ssa^ 
M4 = S4a^. 



(22) 
(23) 



4. Comparison to A^-body Simulations 

4.1 Simulation Methods 

We can check the range of validity of our analytic approximations by comparing their predic- 
tions to the results of fully non-linear, cosmological A''-body simulations. Here we examine simu- 
lations in which the initial conditions are Gaussian with a power law power spectrum, P{k) = /c" 
with n = — 1. We have also analyzed simulations with n = 0; we will discuss these results briefly 
in §4.4. We assume an Einstein-de Sitter {Cl = 1) background cosmology. 

Our simulations use a particle mesh (PM) A^-body code written by Changbom Park. The code 
is described by Park (1990; see also Park & Gott 1991). Its performance has been tested against 
analytic solutions for non-linear pancake collapse and against other PM codes, P^M codes, and 
tree codes. Since the perturbation theory results described above must apply for sufficiently low 
fluctuation amplitudes, comparison of simulations to these results gives us an opportunity to test 
the PM code in a new regime, that of weakly non-linear clustering. 

We wish to examine behavior over a wide dynamic range, so we employ large simulations that 
evolve 100'^ particles on a 200^ force mesh. The code uses a staggered-mesh technique (Melott 1986) 
to achieve higher force resolution (by about a factor of two) than a conventional PM code. We 
advance the particle distributions to expansion factor a = 1/128 via the Zel'dovich approximation, 
then use the PM code to integrate to a = 1 in 127 timesteps of equal Aa. We obtain output at 
expansion factors a = 1/8, 1/4, 1/2 and 1. To analyze the final density fields, we first bin the 
particles onto a 100^ grid using a cloud-in-cell (GIG) weighting scheme, then smooth this field by 
Fourier convolution with Gaussian filters of varying smoothing lengths . We normalize the initial 
power spectrum so that at the final output time (a = 1) the r.m.s. fluctuation predicted by linear 
theory on a Gaussian smoothing scale of = 2 cells is cr = 2. At each output time we analyze 
the density fields with smoothing lengths Vg = 2, 4, and 8 cells. Each of these cells is twice the 
linear size of the cells used for force calculations in the simulations, so the effective gravitational 
softening length is quite a bit smaller than even our smallest smoothing scale. We have run eight 
simulations with independent initial conditions. 

The great advantage of adopting initial conditions with power law power spectra is that one 
can check the reliability of the numerical results by looking for self-similar behavior (see discussion 
by Efstathiou et al. 1988). Neither the form of the initial spectrum nor the $1 = 1 cosmology 
introduces a preferred scale, so at any time only the amplitude of fluctuations is available to define a 
characteristic radius. Statistical properties of the density field smoothed at a particular scale should 
depend only on the r.m.s. fluctuation amplitude on that scale, regardless of whether structure on 
that scale is linear, weakly non-linear, or strongly non-linear. For $1 = 1 and n = — 1, the r.m.s. 
linear fluctuation amplitude is proportional to the expansion factor and inversely proportional to 
the smoothing length, so with our normalization (7(0, r^) = Aa/rg- (If we adopt the standard 
procedure of matching the simulations' mass fluctuations to observed galaxy count fluctuations, 
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then the imphed physical size of the simulation cube is ~ 60a~^ Mpc on a side.) To the extent 
that the simulations are correct, we expect statistical results to depend on a and Vg only through the 
ratio a/rs, and our choices of output times and smoothing lengths provide a number of "degenerate" 
combinations with which to test for this self-similar scaling. Numerical parameters like the force 
resolution, particle density, and box size remain fixed in simulation units, so numerical artifacts 
that reflect the simulation's finite dynamic range should violate this scaling. For the measures that 
we examine in this paper, our simulations obey the self-similar scaling extremely well, except for 
some modest finite- volume effects, which are noticeable at the largest smoothing length, = 8. 

Two features of our simulations are essential to obtaining the excellent agreement between 
A'^-body and pcrturbative results illustrated below. The first is a high initial redshift, so that 
the physical volume represented by the simulation grows by a substantial factor before the first 
output time. If the expansion factor is too small, then the A^-body results will contain transients 
that reflect the use of the Zel'dovich approximation to generate initial conditions. The Zel'dovich 
approximation yields incorrect values for the moments of the density field, primarily because it 
conserves momentum only to first order. For skewness this failure is relatively modest; with a 
top hat filter and P{k) cx fc" power spectrum, for instance, the Zel'dovich approximation yields 
S3 = 4 — (n + 3) instead of the value 6*3 = 34/7 — (n + 3) obtained from rigorous, second-order 
perturbation theory (JBC). However, the approximation deteriorates catastrophically for higher 
order cumulants (Grinstein & Wise 1987; JBC). Furthermore, it fails badly even at the skewness 
level when applied to the divergence of the velocity field; with a top hat filter and power law 
spectrum, the Zel'dovich approximation predicts 830 = {9^)/ {6'^)'^ = Q'^'^{n+1), while perturbation 
theory yields S30 = ri°-^(n — 5/7), i.e. a skewness of the opposite sign for — 1 < n < 5/7 (BJDB). 
We therefore believe that calculations of the PDF based on the Zel'dovich approximation should 
be treated with caution, though they may yield a useful qualitative picture in some cases. 

The second essential feature of these simulations, for our application, is the rather fine (200^) 
force mesh used during dynamical evolution. We first carried out these experiments using a 100^ 
force mesh, and we found ~ 20% discrepancies between the iV-body and analytic results for skewness 
of the density field at low variance, along with comparable violations of self-similar scaling in the 
simulations themselves. The high-resolution mesh is needed precisely because we are investigating 
weakly non-linear clustering. PM codes make significant errors in the linear evolution of Fourier 
modes with wavelengths of a few mesh cells (Bouchet, Adam &; Pellat 1985). However, once 
clustering becomes fully non-linear, structure on the mesh scale is determined mainly by the collapse 
of modes that initially had large wavelengths, and which were therefore evolved accurately through 
the linear regime (Little, Weinberg & Park 1991; Moutarde et al. 1991). We thus expect reliable 
results from a PM code fairly close to the mesh scale provided that the r.m.s. fluctuation exceeds 
unity on the scale of a few mesh cells. Most AT-body studies operate in this regime, and we believe 
that it is this transfer of power from large scales to small scales that accounts for the relatively 
good agreement between different types of cosmological iV-body codes found by Weinberg et al. 
(unpublished comparison). Our present investigation of the weakly non- linear regime requires a 
high-resolution mesh because we are interested in the regime of small fluctuation amplitudes. Note 
that it is specifically the mesh used to compute forces during dynamical evolution that is relevant 
to these considerations, and that a staggered-mesh PM scheme offers significant advantages over a 
traditional PM scheme (Melott 1986; Melott, Weinberg h Gott 1988; Park 1990). 

4.2 Moments and PDF of the Density Field 

Figure 2 compares moments of the TV-body density fields to the predictions of second-order 
perturbation theory. The solid line in the lower panel shows the perturbative relation M3 = [5^) = 
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Figure 2 — Evolution of the asymmetry measures Ms = {5^) and M3 = {5\5\) . Logarithms are base- 10 in this and 
all subsequent figures. In the lower panel, the solid line shows the prediction of second-order perturbation theory, 
= S^cr'^, using the value S3 = 3.47 appropriate to an n = —1 power spectrum and Gaussian smoothing filter. The 
dotted line shows the relation computed from the second-order Edgeworth approximation, M3 = Sau^. Points show 
measurements from the density fields of the AT-body simulations, with smoothing lengths of 2, 4, and 8 cells (circles, 
triangles, and squares, respectively). We use the same set of symbols for M3 and M3, but they can be distinguished 
by their close fits to the corresponding analytic predictions. For closer inspection, the upper panel plots the ratios 
Ma/cr* (top points) and Mz/a^ (bottom points), with horizontal lines representing the analytic predictions. Error 
bars mark the la numerical uncertainty, i.e. the run-to-run dispersion of the eight independent simulations divided 
by 7^/^. Perturbative and iV-body results agree to within this uncertainty when a is small, as expected. iV-body 
results for M3 remain remarkably close to the perturbation theory prediction even when a = 2. There are no free 
parameters to either of these "fits." 

S^(T^, where we have used the value = 3.47 appropriate for an —1 power spectrum and 
Gaussian smoothing filter (JBC). The dotted line shows the relation M3 = {5\5\) = S^a^, derived 
in §3. Points show corresponding results from the A^-body density fields (averaged over the 8 runs), 
with Gaussian smoothing lengths of 2 cells (circles), 4 cells (triangles), and 8 cells (squares). We 
use the same set of symbols for both M3 and M3, but the two quantities can be distinguished 
because of their close fits to the corresponding analytic results. Concentrating first on the A^-body 
results, we see that the circles and triangles agree almost perfectly within their range of overlap, 
indicating that at these smoothing lengths the simulations obey the expected self-similar scaling. 
Values for = 8 (squares) agree quite well, but all three moments have slightly lower amplitudes, 
because the smoothing length is large enough that the absence of power on scales larger than the 
fundamental mode of the simulation cube has become a significant effect. 

The agreement between the AT-body results and the second-order prediction for M3, in a plot 
spanning six orders of magnitude, is rather spectacular. For closer inspection, the top set of points 
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in the upper panel of Figure 2 presents the ratio Ms/o"'^ as a function of cr, and the horizontal solid 
line shows the perturbation theory prediction, M^/a^ = 3.47. Error bars mark the la uncertainty 
in the mean from the eight A'^-body simulations, i.e. the run-to-run dispersion in the value of this 
ratio divided by 7^/^. These error bars represent the uncertainty of our numerical estimates, and 
they increase with smoothing length because the number of independent smoothing volumes per 
simulation decreases. The analytic and numerical results agree to within the statistical error of the 
A?^-body simulations for cr < 1/2. For higher a the N-hoAy results climb above the perturbation 
theory prediction, but they still agree to within 15% at a = 1 and 25% at a = 2, well into the 
regime where second-order perturbation theory should break down. Similar behavior has been 
seen in observational data (Bouchet et al. 1993 and references therein) and in other numerical 
simulations (Lucchin et al. 1993 and references therein). In many cases the near proportionality 
between M3 and extends still further, into the strongly non-linear regime. On the observational 
side, proportionality of moments is implied by the well-known result that the reduced 3-point 
correlation function at small separations can be expressed accurately as a sum of pairwise products 
of the two-point correlation function, (^123 = Q{^i^2 + + CsCi) (Peebles 1980). Among the 
A'^-body results, including our own, there is a rough consensus that the ratio M^/a'^ rises somewhat 
above the perturbation theory value when a exceeds unity, but the change is a modest one. While 
perturbation theory shows that M3 should scale with cr^ when clustering is weak, there is as yet no 
fundamental explanation for the continuation of this relation into the strongly non-linear regime. 

^ Figure 2 also shows agreement between the A?^-body results and the Edgeworth calculation of 
M3 for small a. However, in this case the analytic predictions begin to separate from thenumerical 
results more quickly. The pcrturbative calculation overestimates the A^-body values of M3 by 15% 
at C7 = 1/2, 30% at cr = 1, and 60% at a = 2. Even these errors are not so large when one recalls 
that the predicted value of M3 is growing by a factor of 8 for each factor of 2 increase in cr. Positive 
<md negative fluctuations grow at different rates in the non-linear regime, and since the M3 and 
M3 moments weight extreme values differently, it would be virtually impossible for perturbative 
calculations of both quantities to remain accurate once cr approaches or exceeds one. 

In Figure 3 we plot the width measures M4 and M4 against cr. The solid line in the lower 
panel shows the prediction of third-order perturbation theory for the fourth-order cumulant: M4 = 
((5^) — 3cr^ = S'4cr^. Points show the Ai'-body results for Ts = 2 (circles) and = 4 (triangles). The 
fourth moment is very noisy at = 8, and we do not show results for this smoothing length. In 
the upper panel, the top set of points shows the ratio of M4 to a^, with error bars representing 
the Icr uncertainty from the simulations and the horizontal solid line representing the prediction 
M/^/a^ = constant, from perturbation theory. We do not have an analytically calculated value of ^4 
for this spectrum and smoothing filter, so we have treated it as a free parameter and chosen a value 
5*4 = 20 that provides a good eye-fit to the low-cr points in Figure 3. This value is in reasonable 
accord with GGRW's computation for a CDM power spectrum; they find S'4 20 on a smoothing 
scale where the slope of the CDM power spectrum is n ~ — 1. 

The dotted line in the lower panel of Figure 3 displays the relation — M4 = (2/7r)^/^cr— {\5\) = 
—S^a^ predicted by the Edgeworth approximation (equation 23; note that ^4 is negative). We 
compute the value of 54 using S'4 = 20, determined from the fit to the M4 vs. plot, and 
S3 = 3.47, determined analytically, so there is no freedom to adjust the height or slope of the 
dotted line. Points show the corresponding TV-body results. In the upper panel, the bottom 
set of points and the dotted line show the ratio — M4/cr'^ for the simulations and the Edgeworth 
approximation, respectively. We have multiplied all the values and the error bars by a factor of 30 
in order to make them clearly visible on this plot. The results of Figure 3 are similar to those of 
Figure 2. For both M4 and M4, the A?^-body points follow the power laws predicted by perturbation 
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Figure 3 — Evolution of the width measures M4 = (5^) — 3(J^ and — M4 = (2/7r)^/^(T — In the lower panel, the 

solid line shows the prediction of third-order perturbation theory, M3 = 84^(7^'. The dotted line shows the relation 
computed from the third-order Edgeworth approximation, — M4 = —S4a^ (note that 6^ is negative). Points show 
measurements from the density fields of the AT-body simulations, with smoothing lengths of 2 cells (circles) and 4 
cells (triangles). The upper panel displays the ratios M^ja^ and 30 x ( — M4/(t^), with the factor 30 chosen to make 
the results easily visible on this plot. We do not have an analytically calculated value of 54, so we have chosen a 
value 54 = 20 that provides a good eye-fit to the A/"-body results for M4 at low cr. This constant is the only free 
parameter in these "fits" ; the slopes of the power laws are determined by perturbation theory, and the value of S4 
is fixed by the choice of S4. 



theory when a is small. The cumulant M4 remains close to the perturbation theory prediction even 
for cr ~ 1 — 2. The Edgeworth approximation for M4 stays fairly accurate up to cr = 1, but it 
overestimates the numerical results significantly at cr = 2. 

Figure 4 brings us to the central issue of this paper, the overall shape of the PDF. The solid 
lines show PDFs of the iV-body density fields for three different values of the r.m.s. fluctuation 
amplitude, (7 = 1/8 (top panels), cr = 1/4 (middle panels), and cr = 1/2 (bottom panels). Left hand 
panels plot p(z^) against v = 5/a and emphasize behavior near the peaks of the distributions. Right 
hand panels plot log^o^^ against u in order to display the tails of the distributions. We average the 
results from the eight A^-body density fields analyzed with smoothing length = 4; the curves for 
different values of a are obtained from output expansion factors a = a. We obtain identical results 
from other combinations of a and Vg that have the same a, except for minor finite-volume effects 
on the extreme tails of the distributions. 

The dotted lines in Figure 4 show the second-order Edgeworth approximation to the PDF 
(equation 13), i.e. including only the Gaussian and skewness terms of the Edgeworth series. The 
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Figure 4 — Probability distribution functions (PDFs) evolving from Gaussian initial conditions. Left hand panels 

are linear plots, emphasizing behavior near the peak of the PDF. Right hand panels are logarithmic, emphasizing 
the tails. Solid lines show PDFs measured from the smoothed density fields of the A'^-body simulations when the 
r.m.s. fluctuation is 1/8 (top), 1/4 (middle), and 1/2 (bottom). Dotted lines show the second-order Edgeworth 
approximation to the PDF, equation (13). Dashed lines show the third-order Edgeworth approximation, equation 
(14). Both approximations work well for (7 = 1/4 and begin to break down when a = 1/2. 

approximation does very well at a = 1/8, and it remains accurate at cr = 1/4 except that the 
positive tail falls too rapidly for > 4. At a = 1/2 the approximation is breaking down, though 
it still oscillates around the true PDF in a reasonable way. The dashed lines show the effect of 
including the third-order term of the Edgeworth expansion (equation 14). We use the value 5*4 = 20 
estimated from the moment ratios in Figure 3. This third-order approximation is somewhat more 
accurate than the second-order approximation at each value of a. This is the behavior that we 
expect, since our approximation is a power series expansion in cr. The third-order term provides 
higher accuracy, and it slightly extends the useful range of the Edgeworth expansion. However, the 
third-order approximation begins to break down fairly soon after the second-order approximation 
fails; the Edgeworth series is a powerful approach for describing the first deviations from Gaussian 
fluctuations, but it cannot take one into the deeply non-linear regime. Also, as Figure 4 illustrates, 
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Figure 5 — Evolution of — M35) = — {0^} and —M^g = —{6\8\), where 6 = V ■ v/Hq is the divergence of the pecuhar 
velocity field. The format is the same as Figure 2. We use the value 5351 = —2.19 appropriate to an n = —1 power 
spectrum and Gaussian smoothing filter, so there are no free parameters to either of these "fits." 



the Edgcworth approximation to tlie PDF is not positive-definite, though it does integrate to unity 
at any order, and it is positive-definite when a <C 1. We have not investigated the effects of including 
higher-order terms of the Edgeworth series. If one uses the Gram-Charher series (1) instead of the 
Edgeworth series (12), then adding the third-order term degrades the fit to the A'^-body PDF. This 
failure is not surprising, since the Gram-Charlier series is not a proper asymptotic expansion. 



4.3 Moments and PDF of the Velocity Divergence 

We now turn our attention to the velocity field, or, more specifically, to its divergence. Per- 
turbative calculations describe a smoothed velocity field in which each volume element has equal 
weight. However, AT-body simulations have a finite number of particles, not a continuum, and the 
simulation velocity field is known only at the particle locations. When the clustered particle distri- 
bution is binned onto a grid, some cells may be empty, but this docs not mean that the velocities 
in those cells are zero, just that they are undetermined. It is straightforward to define a mass- 
weighted, smoothed velocity field as the ratio of the smoothed momentum field to the smoothed 
density field, but the results may not agree with perturbative calculations because of the difference 
in definitions. [A velocity field defined on a grid by averaging the velocities of the particles in each 
cell is, in fact, a mass-weighted velocity field with a cubic cell smoothing kernel.] To get around 
this problem, we define the velocity field smoothed on scale by first computing a mass-weighted 
velocity field on a grid with a Gaussian smoothing length rs/2, then smoothing this field in a 
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Figure 6 — Evolution of M^^g and — Af4g, in the same format as Figure 3. In tlie upper panel, we multiply the 
ratios —Mi^g/a'g (bottom set of points and horizontal dotted line) by a factor of 20, to make the results easily visible 
on this plot. We do not have an analytically calculated value of 54^, so we have chosen a value S^g = 6.0 that 
provides a reasonable ^e-fit to the iV-body results. However, the results alone suggest a somewhat lower value 
{Si w 4.5), while the M^g results alone suggest a somewhat higher value {Si « 6.5). The slopes of the power laws 
are determined by perturbation theory, and the relation between S^g and S^g is determined by equation (21), so the 
value of Sig is the only free parameter in these "fits." 

volume- weighted way with a smoothing length r'^ = (r^ — The first step assigns sensible, 

non-zero velocity values to all cells, but since smoothing lengths add in quadrature, the second, 
volume- weighted smoothing dominates the final result, and the combined smoothing length is equal 
to Vg- When density fluctuations are small, this procedure is equivalent to simple, volume- weighted 
smoothing of the velocity field. Once the smoothed velocities have been calculated, we compute 
the PDF and moments of the divergence field 9 = V -v/Hq. In the linear regime, —9 is equal to the 
density contrast 5, but the non-linear evolution of the two fields is different even at second order. 

Figure 5 plots —M^e = —{9^) and — M30 = —{9\9\) against ae = (6*^)^/^, in the same format 
as Figure 2. The PDF of the velocity divergence has negative skewness, because in the non-linear 
regime the inflows to high density regions are faster than the outflows from low density regions. For 
the analytic predictions we use the coefficient S30 = —2.19 appropriate to the velocity divergence 
field (BJDB). Once again, the A^-body results obey the expected self-similar scaling except for 
minor finite-volume discrepancies at = 8. The horizontal spacing between points decreases 
steadily, indicating that ag is growing more slowly than predicted by linear theory. By the final 
time, the r.m.s. fiuctuations at 2 and 4 cells are ae{2) = 1.0 and ere (4) = 0.68, compared to the 
linear theory values of 2.0 and 1.0. The r.m.s. fluctuation of the density fleld, on the other hand, 
grows at almost exactly the linear rate (see Figure 2). 
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Figure 7 — Evolving PDFs of the velocity divergence, in the same format as Figure 4. Solid lines show PDFs 

measured from smoothed velocity fields of the A^-body simulations, when the r.m.s. fiuctuation of the velocity 
divergence is 0.24 (top), 0.43 (middle), and 0.68 (bottom). Dotted and dashed lines show the second- and third- 
order Edgeworth approximations, respectively. 



Figure 5 demonstrates agreement between the iV-body and perturbative calculations for low 
(T0. With the n = — 1 power spectrum and Gaussian filter used here, the Zel'dovich approximation 
predicts a skewness that is lower than these values by nearly a factor of 8 (BJDB). The second- 
order predictions remain close to the A^-body results even as ag approaches 1, overestimating the 
numerical values of -M30 and -Mg^ by about 15%, 20%, and 25-30% for ae = 0.42, 0.68, and 
1.0 respectively (corresponding to linear theory values of ae = 0.5, 1.0, and 2.0). Figure 6 shows 
corresponding results for M^e = {0^} — 3crg and —M40 = {2/Tr)^/^ae — {\9\), in the same format 
as Figure 3. In the upper panel, we multiply the ratios —M^g/ag by a factor of 20 to make them 
clearly visible on this plot. We do not have an analytically calculated value of 84^0, and we have 
chosen the value S^q = 6 on the basis of an eye-fit to the A^'-body points in Figure 6. Small changes 
in the value of S40 induce large logarithmic changes in the value of S40 oc S^q — S40 (equation 21), 
because it is in a critical region near zero. If we were fitting the M40/a^ ratio alone, we would 
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probably choose 84$ ~ 6.5, but the results for seem to indicate a lower value, £'4^1 ~ 4.5. We 
do not know the cause of this mild discrepancy, but we would guess that it reflects inaccuracies and 
statistical uncertainties in the numerical results, the significance of which is amplified by being in 
the critical region 84$ ~ Sfg. An analytically computed value of S40 might clarify this issue, but 
the required integrations are rather forbidding. 

Figure 7 compares the Edgeworth approximation for the PDF of 9 to the PDFs measured 

from the iV-body simulations, in a format similar to Figure 4. Solid lines show the A^-body PDFs 
with a smoothing length of 4 cells at o = 1/4, 1/2, and 1 (top to bottom). The values of the 
r.m.s. fluctuation ae are listed in each panel. Dotted lines shows the second-order Edgeworth 
approximation, equation (13) with 83$ = —2.19. Dashed lines show the third-order approximation, 
equation (14) with S^g = 6. The PDF of the velocity divergence develops a non-Gaussian shape 
more slowly than the PDF of the density fleld, and the second-order Edgeworth approximation 
remains accurate further into the non-linear regime, as one can see by comparing to Figure 4 
(where the r.m.s. fluctuations are actually lower). It comes as no great surprise that the accuracy 
of the second-order Edgeworth expression is closely related to the magnitude of the skewness, jS'sO"!, 
or that the approximation begins to break down when |53C7| exceeds one. The third-order term of 
the expansion generally improves the behavior of the negative 9 tail, but overall this term appears 
to be less useful for the velocity divergence than it is for the density field (compare to Figure 4). 

4.4 Spectral Dependence of the PDF 

We have carried out all of the above comparisons for similar simulations with a white noise 
(n = 0) initial power spectrum. We do not show the results here, but the agreement between the 
perturbative and A''-body calculations is at least as good, and in some cases even better, provided 
that one uses the values 5*3 = 3.14 and 8se = —1.67 appropriate to an n = spectrum (JBC; 
BJDB). This leads us to an interesting theoretical point. Kofnian et al. (1993) compute the density 
PDF that evolves from Gaussian initial conditions using the Zel'dovich approximation. In their 
calculation the shape of the PDF depends only on the r.m.s. fluctuation amplitude, a, and it is 
independent of the shape of the power spectrum. Strictly speaking, Kofman et a/.'s technique 
describes the PDF of an unsmoothed density field evolving from smoothed initial conditions, and 
if we used the second-order Edgeworth approximation for this PDF we would also find a spectrum- 
independent result, since the coefficient 5*3 = 34/7 is independent of the power spectrum in the 
absence of smoothing. For the case of a smoothed final field, the second-order Edgeworth PDF 
depends on both the r.m.s. fluctuation and the shape of the power spectrum, but only through 
the combination 8sa. This single parameter tells us how to relate the predicted PDFs for different 
power spectra, and how to relate the PDF of the density field to that of the velocity divergence. 

We can conjecture that the shape of the PDF may continue to depend primarily on the param- 
eter 83a even when the second-order Edgeworth approximation begins to break down. In abstract 
form, we can write this conjecture as 

p[iy;n;a{rs)]=p[iy;83{n)a{rs)] ; (24) 

i.e. the shape of p(z^), which could in principle depend on the spectral slope n and on the fluctuation 
amplitude a at the smoothing scale r^, in fact depends on these parameters only through the 
combination 83{n)a{rs)- In the same notation, we can express the scaling proposed by Kofman 
et al. (1993) as 

p[u; n; a{rs)] = p[iy; o-(rs)]. (25) 
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Figure 8 — Comparison of density PDFs for different initial power spectra. The heavy solid lines show the PDF 
measured from iV-body simulations with a white noise (n = 0) initial power spectrum. The r.m.s. fluctuation on the 
smoothing scale is u = 0.80, so S-^a = 3.14 x 0.80 w 2.5. The dotted lines (often obscured by the heavy solid lines) 
show the PDF of the n = —1 simulations at a scale where a = 0.73 and Sacr = 3.47 x 0.73 ~ 2.5. Even though the 
second-order Edgeworth approximation (indicated by the light solid lines) has failed rather badly, the PDFs of the 
n = and n = — 1 models are a nearly perfect match on the scales where S^cr is the same. The dashed lines show 
the PDF measured from the n = —1 iV-body simulations at a scale where a = 0.80 and Sz<J = 3.47 X 0.80 « 2.8. 
The shape of the evolved PDF is more nearly a universal function of Sza than it is of u alone. 

Figure 8 presents a test of the conjecture (24). The heavy sohd hues in each panel (hnear on the left, 
logarithmic on the right) show the PDF measured from our n = simulations at a = 1/2, with a 
smoothing length of 2 cells. The r.m.s. fluctuation on this scale is cr = 0.80. The dotted lines, which 
are mostly obscured by the heavy solid lines, show the PDF measured from our n = — 1 simulations 
at a = 1/2, with a smoothing length of 2.82 cells. The corresponding a is 0.73. The values of S^a 
for these two sets of density fields are very nearly equal: 3.14 x 0.80 3.47 x 0.73 2.5. The light 
solid curve shows the second-order Edgeworth approximation for this value of Sscr. We see that 
the PDFs of the two fields match very closely, even though neither of them agrees well with the 
analytic approximation. The dashed line shows the PDF of the a = 1/2, n = —1 simulations at 
a smoothing length of 2.55 cells, the scale where the r.m.s. fluctuation a = 0.80 matches that of 
the n = density fields. While this PDF is reasonably close to that of the n = simulations, it 
is clear that the shape of the PDF depends on the spectral slope as well as on the value of a, and 
that equation (24) offers a more accurate description of the PDF than equation (25). Analytically, 
we can see that equation (24) will hold at the level of the third-order Edgeworth approximation 
(equation 14) if and only if the ratio 5*3 / 5*4 is independent of n. Bernardeau (in preparation) shows 
that this ratio is indeed nearly independent of n for a top hat smoothing filter. 

5. Biased Galaxy Formation 

Perturbation theory describes the evolution of the mass distribution, but observations probe 
the distribution of galaxies. There are both observational and theoretical reasons for thinking that 
galaxies do not evenly trace the large-scale mass distribution. On the observational side, we know 
that elliptical and spiral galaxies have different clustering properties; it is clear that the two galaxy 
types cannot both trace the mass independently, and there is no particular reason to expect that the 
union of the two classes does trace the mass. Theoretically, we know that the collapse epoch of a 
galaxy scale perturbation will depend on the background density, so the history of perturbations will 
vary with environment, and the efficiency of galaxy formation may vary correspondingly. Numerical 
simulations that include gas dynamics indicate that galaxy formation is at least somewhat biased 
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towards regions of high background density (Cen & Ostriker 1992; Katz, Hernquist & Weinberg 

1992) . 

For the mass distribution, second-order perturbation theory tells us that {5^) = S^a"^, if the 
primordial fluctuations are Gaussian. Is there a corresponding relation for the galaxy density 
contrast Sg7 If we adopt the simplest mathematical relation between the galaxy and mass density 
contrasts, the linear bias model dg = bS, then the answer is obvious: {6g) = S^gcr^, where Szg = 
S:i,/h. However, there is no theoretical motivation for the linear bias model, and while it may be a 
useful approximation for some purposes, it seems risky to assume a symmetric bias relation in order 
to compute a measure of asymmetry in the probability distribution. Indeed, one might worry that 
allowing a non-linear relation between galaxy density and mass density would destroy the simple 
relation between [5^) and predicted by perturbation theory, but this is not the case, as we shall 
now see. 

Instead of a linear bias, let us adopt the much looser assumption that the smoothed galaxy 
density is a non-linear but local function of the smoothed mass density (see discussion by Coles 

1993) , 

5g = m. (26) 

The second-order Taylor expansion for Sg is 



32 c"^, wJ 

6 = /'(0), h = nO), a^ = {S'). (27b) 



dg = bS + ^ib"^ - -^hcr^, where (27a) 



The last term on the right-hand side of equation (27a) is required to ensure that {5g) = 0. It is 
straightforward to use the expansion (27) to compute and (Sg) to O(cr^), making the substitutions 
{6^) = Sscr"^ and {6"^) = 3a^ + S4a^ where appropriate. The result is 

= 6V^ Msg = {S^ = {Ssb^ + 3bH2)(T\ (28) 

making the moment ratio for the galaxy distribution 

_Msg_Ss 3b2 

For linear bias, 62 = 0, we recover the earlier result Ssg = Ss/b. However, the value of Ssg depends 

sensitively on the shape of the bias function through the second-derivative term 362/^^- This shape 
dependence may explain why the value of Ssg measured by Bouchet et al. (1993) from the IRAS 
redshift survey is rather low, S^g ~ 1.5. IRAS galaxies are underrepresented in rich clusters relative 
to optical galaxies, so the "bias function" that applies to IRAS galaxies may have negative curvature 
(negative 62)) pushing S^g below the value of 5*3 for the mass. Evidence that cluster-avoidance is 
an important effect comes directly from Bouchet et a/.'s (1993) analysis; they find that double- 
counting the IRAS galaxies in rich cluster cores, a small fraction of the total sample, more than 
doubles the measured value of S^g- 

The most important implication of equation (28) is that {5^) will be proportional to on 
scales where ag is small, provided that the linear mass fluctuations are Gaussian and that the 
galaxy density is a local function of mass density. Fry & Gaztaiiaga (1993) have independently 
derived equation (28), and they have generalized the result in an important way: by expanding 
5g = f{S) in higher-order Taylor series, they show that a local biasing function preserves all of the 
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Figure 9 — Evolution of M^g = {5g) and M^g = {5g\Sg\), where 5g is the contrast of the galaxy density field, related 
to the mass density field by the "biasing function" (30). Format is the same as Figure 2. The solid and dotted lines 
show the analytic relations Msg = SsgCfg and Mzg = Ssgag. Circles, triangles, and squares show results from the 
biased A^-body density fields at smoothing lengths of 2, 4, and 8 cells, respectively. The value of Szg is computed 
from equation (29), so there are no free parameters to either of these "fits." When ag is small, biasing preserves the 
form of the moment relations predicted by perturbation theory, and with this biasing function, which is derived from 
hydrodynamic cosmological simulations, the perturbative relation M^g = Szg(Tg remains accurate even at ag « 3.5. 



moment relations predicted by perturbation theory (equation 11), in the hmit of smaU fluctuation 
amplitude. One can therefore test the hypotheses of Gaussian primordial fluctuations and local 
biasing by examining the moments of the galaxy count distribution on large scales. 

We can illustrate these analytic arguments and get a sense of how they extend into the non- 
linear regime by considering the bias function proposed by Ccn & Ostrikcr (1993; hereafter CO), who 
incorporate a simple but physically motivated recipe for galaxy formation into their hydrodynamic 
simulations of the cold dark matter model. CO fit the relation between galaxy density and mass 
density in their simulations to a non-linear functional form, 

log(l + 6g) = A + Slog(l + 5) + C[log(l + 5)f. (30) 

We can apply this transformation directly to the smoothed mass density fields of our AT-body 
simulations and compute the resulting moments. In their simulations, CO find that the constants 
B and C depend weakly on smoothing scale. We take the values B = 1.5 and C = —0.14 that CO 
find for a Gaussian filter scale of 5^"-^ Mpc; at any output time and smoothing scale, the constant 
A is then determined by the requirement that {5g) = 0. The details of this biasing procedure 
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Figure 10 — Evolution of M^g and —M4g, in the same format as Figure 3. Solid and dotted lines in the lower 
panel show the analytic predictions M4g = S4gag and —M^g = —SAga^. Circles and triangles show results from 
the biased iV-body density fields at smoothing lengths of 2 and 4 cells, respectively. The value of S^g is computed 
analytically using the value of ^4 measured for the sjrnulation mass density fields, so there are no free parameters 
to these "fits." In the upper panel, the values of —M^g/a^ (bottom set of points and horizontal dotted line) are 
multiplied by a factor of 30, to make them easily visible on this plot. 



should not be taken too seriously, in part because the resolution of the CO simulations themselves 
is only ^ 0.5/i~-^ Mpc, which is rather low for inferring rates of galaxy formation. Nonetheless, this 
biasing scheme reflects the sort of relation between galaxy and mass density that one might expect 
in a rather broad class of theoretical scenarios. (The simulations of Katz et al. [1992] and Evrard, 
Summers &; Davis [1993] have much higher spatial resolution in galaxy forming regions, but in each 
case the simulation cube is ~ 10/i~^ Mpc, too small a volume in which to measure a meaningful 
biasing function.) 

Figure 9 plots M^g and M^^g against Gg for our A''-body simulations, where 5g and 5 are related 
by equation (30). The format is the same as Figure 2. The solid and dotted lines in the lower panel 
show the perturbation theory predictions, M^g = SsgfJg and M^g = Ssgffg, respectively. The upper 

panel plots the ratios M^g/Ug and Msg/ag for the simulations and the analytic approximations. 
We compute Ssg = 3.15 from equations (29) and (30) in the limit of (5^) 0; in this limit, the 
constant A in the biasing function (30) goes to zero. The results are strongly reminiscent of those 
for the mass distribution, shown in Figure 2. When cjg is small, the A^-body points sit precisely 
on the perturbation theory lines, as expected. As ag approaches unity, the analytic approximation 
for M^g begins to fail, but M^g follows the perturbation theory prediction remarkably closely up to 
the last point computed, where ag ^ 3.5. This continuing agreement is not guaranteed by equation 
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(28), since the derivation of that equation is based on a second-order Taylor expansion which is no 
longer valid when ag <^ 1. 

Figure 10 plots M^g and against ag, in the same format as Figure 3. The solid line in the 
lower panel shows the perturbative relation = S^gag. We compute S^g = 16.9 using equation 
(10) of Fry &; Gaztaiiaga (1993), assuming 5*4 = 20 for the mass (from Figure 3). The dotted line 
shows the Edgeworth approximation M^g = Siga^, with 5*43 computed from equation (21). The 

upper panel displays the ratios M4g/ag and Mig/a^, with the latter multiplied by a factor of 30 for 
visibility. The analytic and numerical results agree when ag is low, as they should. The Edgeworth 
approximation of M4C, overestimates the A^-body value at larger ag, but the value of M^g stays 
close to the relation predicted by perturbation theory up to the last point, ag « 3.5. We have 
also compared the full PDFs of the biased density fields to those computed from the Edgeworth 
expansion. The results are similar to those shown in Figure 4 for the mass distribution: good 
agreement for ag ^ 1/2, and poor agreement beyond. 

If the primordial mass fluctuations are non-Gaussian, e.g. if they have intrinsic skcwness or 
kurtosis, then the linear term of the biasing function will transfer these intrinsic moments to the 
galaxy fluctuations, and moments of the galaxy counts will not obey the hierarchical relations of 
equation (11). Galaxy counts may also violate equation (11) if the galaxy density is not tightly 
coupled to the local mass density. One could imagine, for instance, that the galaxy density obeys 
Sg = f{5) on average but with substantial scatter about the mean trend. Scatter might arise if the 
efficiency of galaxy formation is sensitive to the pressure of the local intergalactic medium or to ion- 
izing radiation from nearby quasars. In the limit where scatter about the mean relation overwhelms 
the trend predicted by the mean relation itself, it is clear that moments of the galaxy distribution 
will reflect the moments of the "scatter function" rather than moments of the underlying mass 
distribution, and there is no reason to expect these moments to obey the special relations implied 
by equation (11). For example, Frieman & Gaztaiiaga (1993) show that the "cooperative galaxy 
formation" scheme, proposed by Bower et al. (1992) as a possible way to reconcile the standard 
CDM model with the galaxy angular correlation function of the APM survey (Maddox et al. 1990), 
predicts a strong shift in the relation between skewness and variance of galaxy counts at scales 
~ 10 — 20/i~^ Mpc. Precise observational confirmation of hierarchical relations between moments 
of galaxy counts would provide evidence in favor of Gaussian primordial fluctuations and important 
constraints on the process of galaxy formation itself. 

6. Discussion 

The comparisons in Section 4 provide encouraging news both for the perturbative analytic 
approach and for the A^-body simulations themselves. Most previous tests of cosmological A^-body 
methods have examined either the linear growth of fluctuations or strongly non-linear problems 
like pancake collapse. When the variance is small, the non-linear effects discussed in this paper are 
quite subtle, as evidenced by our plots of dimensionless quantities extending to ~ 10~^ and even 
below. Nonetheless, the moments of our A^-body density fields match the perturbative calculations 
perfectly when the variance is small. These precise measures test the A^-body method in a new 
regime, that of weakly non-linear clustering, and the agreement with analytic theory strengthens 
our faith in the reliability of the simulations. Most cosmological A?^-body studies use a cubic lattice, 
the Zel'dovich approximation, and periodic boundaries to set up initial conditions, just as we 
do. The match to second- and third-order perturbation theory in the weakly non-linear regime is 
significant because it shows that such initial conditions allow A/"-body simulations to settle into the 
correct non-linear solution for the evolution of the mass distribution. Success is not guaranteed by 
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the initial conditions algorithm itself, since the Zel'dovich approximation does not yield the correct 
relations between moments of the density field (Grinstein &; Wise 1987; JBC; BJDB; see discussion 
in §4.1). 

The A^-body simulations confirm the correctness of the analytic calculations, including the 
moment calculations of JBC and BJDB and our new results for PDFs, M3, and M4 of the smoothed 
density field and the smoothed velocity divergence. The simulations show that these results continue 
to hold on large scales even when small-scale clustering is strongly non-linear. This conclusion is 
unsurprising; more remarkable is the fact that the perturbative approximations do not break down 
rapidly in the non-linear regime. In particular, the skewness and kurtosis of the density field stay 
impressively close to the predictions of perturbation theory even when a = 2. The second-order 
Edgeworth approximation to the PDF of the density or velocity divergence remains accurate Tintil 
IS'sO'l, the magnitude of the skewness, reaches one, and the third-order approximation remains 
accurate slightly longer. Calculations of M3 and M^g based on the Edgeworth series match the 
iV-body results to ~ 15% when a ~ 1/2 and ~ 30% when cj ~ 1. We have applied our techniques 
specifically to the case of Gaussian initial conditions. We believe that a similar treatment is possible 
for non-Gaussian initial conditions, though the method requires some modification because a low- 
order Edgeworth expansion may not provide a good description of the linear theory PDF in a 
non-Gaussian model. 

One of the encouraging results of this paper (derived independently by Pry Sz Gaztafiaga 1993) 
is that "biased" galaxy formation preserves the relation between the skewness and variance of the 
density field predicted by perturbation theory, provided that the galaxy density is a local function 
of the mass density. This fact can be demonstrated analytically in the limit of small fiuctuations, 
and once again our tests on numerical simulations show that the predictions of perturbation theory 
continue to hold remarkably well in the fully non-linear regime, at least for the biasing relation 
proposed by Cen &; Ostriker (1993). By studying galaxy counts on large scales, one can learn about 
both the nature of primordial fluctuations and the physics of galaxy formation. 

We have limited the analysis in this paper to the case of = 1. Perturbation theory predicts 
that the S3 coefficient for the density field should have only a very weak dependence on J7 (BJCP). 
However, for the velocity divergence there is a fairly strong dependence, roughly S30 oc 
(BJDB), so predictions for the moment relations (Figures 5 and 6) and the PDF (Figure 7) depend 
significantly on Q. If perturbation theory works equally well for low-fi models, and we have every 
reason to think that it will, then the moments and the PDF of the velocity divergence can be used 
to constrain the density parameter provided that (a) one adopts the hypotheses of gravitational 
instability and Gaussian fluctuations, and (b) one can obtain a reliable estimate of the velocity 
field over a sufficiently large volume. Since this technique does not use the galaxy density field, it 
is independent of biased galaxy formation so long as galaxies provide fair tracers of the large-scale 
velocity field. We address these ideas more fully elsewhere (BJDB). 

A close relative of the velocity divergence technique mentioned above is the "reconstruction" 
method of Nusser & Dekel (1993, hereafter ND), who attempt to recover the PDF of the initial 
density fluctuations from the divergence of the present day velocity field by applying the Zel'dovich 
approximation. They find that the velocity field inferred by POTENT (Bertschinger Sz Dekel 
1989; Bertschinger et al. 1990) is consistent with Gaussian initial conditions if = 1, but not 
if J7 = 0.3. We have two cautionary remarks to make about this conclusion (in addition to the 
caveats listed by ND themselves). First, as discussed in §4.1 and §4.3, the Zel'dovich approximation 
makes large errors in predicting the skewness of the velocity divergence, which is the simplest 
measure of asymmetry in its PDF. Second, the residuals between ND's recovered initial PDF 
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and a true Gaussian, plotted in figure 5 of ND, bear a remarkable resemblance in shape to the 
Gaussian derivatives (p^^^^ and (f)^^'' (see Figure 1). Since these derivatives appear in the Edgeworth 
approximation to the evolved PDF, their appearance in ND's residuals may indicate a systematic 
dynamical failure of their reconstruction method. However, the magnitude of such a failure is 
constrained by the method's reasonable success when applied to A'^-body simulations with known 
initial conditions (ND; Gramman, Weinberg &; Nusser, in preparation). 

With the rapid growth in galaxy redshift and peculiar velocity data, the regime of weakly 
non-linear clustering is becoming increasingly accessible to observations. By comparing PDFs and 
moments of the density and velocity divergence fields to the predictions discussed here, we can 
test the hypothesis that structure in the universe formed by gravitational instability from Gaussian 
primordial fluctuations. In the analysis of galaxy density fields one must introduce assumptions 
about the relation between galaxies and mass, and in the analysis of velocity fields one must 
introduce assumptions about J7, but by examining structure over a variety of scales one can check 
all of the input assumptions for internal consistency. Application of these methods to high quality, 
large volume data sets should therefore teach us a great deal about the formation of galaxies and 
the origin of large-scale structure. 
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Appendix 



Our purpose here is to derive {S\S\), using gravitational instability theory, and the assumption 
that Si is a Gaussian random field. Since 62 = 0{Sl), we can use the expansion 



\S\ = ^{61 +62 + ...y = \6i\ (1 + 62 /5i) + 0{5l) {Al) 

This gives 

5\5\=5i\5i\ + 2 52\5i\+0{5t) . {A2) 

By symmetry, the mean value of the first term above is zero, and the lowest order contribution to 
{5\5\) comes from the second term. To calculate 

(,5|,5|)=2(,52|<5i|), (A3) 

we need the joint probability distribution for 5i and ^2. For density fluctuations in an expanding 
universe, filled with a pressureless non-relativistic fluid, the relevant perturbative solutions are 
given by (Juszkiewicz k, Bouchet 1992) 

(5i = D{t)e{^) , {AA) 

62 = D^{t) + K)e^ + Ve • V$ + (I - Ac)T„^rc,/3] , {A5) 

where we use the Einstein summation convention, t is the cosmological time, x = {x^}, a = 1, 2, 3, 
are comoving coordinates, and D{t) is the linear order fluctuation growth rate (cf. Peebles 1980). 
Without loss of generality, in our remaining calculations below we will always assume D{t) = 1. 
The quantities $ and Tq,i3 arc proportional to the linear order Newtonian potential and the shear 
tensor, respectively: 

= £(x) ; Taf3 = (|<5a/3V2 - V«V;3)$(x) . {A6) 

Here 6^/3 is the Kronecker delta. The parameter k is a slowly varying function of cosmological 
time. For densities in the range 0.05 < < 3, it is well approximated by K[f2(i)] ~ (3/14) 
(BJCP). It is convenient to introduce the Fourier transform 

The power spectrum is defined by the relation 

(£k£k') = i^(/c)<^i?(k + k') , (AS) 

where 6d is the Dirac delta function. The Fourier transform of $ is $k = —^k./k'^- The transforms 
of the potential and density gradients, V$ and Ve , are ikek/^^ and — ikek , respectively. 
The shear tensor can be represented by the Fourier integral 



^«/3(x) = / 77^^ (5<5a/3 - kakp) e"^*''^ , {A9) 



(27r)3/2 

where ka = ka/k. To calculate {5\5\), we need the joint probability density, p(£:, Ve, V$, Tq,^). 
Since e is Gaussian, ^(e, . . .) is a multivariate Gaussian distribution, entirely determined by its 
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covariance matrix. There is no need to calculate this horrifying 12 x 12 matrix explicitly because 
£ is statistically independent from the remaining variables. Indeed, e is correlated only with itself: 

(e^(x)) = / ^(^.^W e-^-^-'y- = J ^ P{k) = , (^10) 

while 

^ikP{k)=o, (-v$)=y ^^p(fc)=o, (All) 

and 

{£T^p)=0. (A12) 
To derive the last expression above, we used equations (A9) and (A7), as well as the identity 

P{k) k^k0^ = , (A13) 

which can be verified easily by integrating over the appropriate angles. Equations (All) and (A12) 
imply 

p{e, V£, V$, r«/3) = p{e) p{Ve, V$, r„^) . {AU) 
Now, let us represent S2 in equation (A5) as a sum of two terms, 

/i = I (1 + k) e^ and Is = Ve ■ V$ + (^ - k) r^p r^p . {Alh) 

The factorization of the joint probability density (eq.[Al4]) allows us to write 

{\e\h) = {\e\){h). (A16) 

The second and higher order corrections to 5i must preserve {5) = because mass is conserved (cf. 
§18 in Peebles 1980). This implies {82) = {h + h) = 0, and therefore 

= -(/,) = -I (1 + K)a^ (AU) 

Now the calculation of {5\5\) is reduced to 

(<5|<5|) = 2 (7i + h)) = 2 h) - 2 {\s\) ih) , {A18) 



and all averages above involve only the marginal distribution p{e) = exp{—e^/2a^) /Vzttct^. The 
contribution from Ii amounts to 



{\e\) = ^= I ee-^'^/^'''' de= J-a , {A20) 



The mean value of |£| equals 

(T\/27r Jo V TT 

so the contribution from I2 is 

'32 



2{\s\){h) = -^^{l + K)a' . {A21) 
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Combining this with eq.(A19), we obtain 



/ QO / 2 





{A22) 



where we use the expression 



^3 =4(1 + ^^) , 



(^23) 



obtained for the unsmoothed field by BJCP. Our above result is in agreement with equation (17) 
in the main text. However, equation (17) is more general than our present result, since its validity 
is not restricted to the unsmoothed field and its particular value of S3. It is valid for arbitrary 
filters and arbitrary power spectra. An attempt to include smoothing would make the calculation 
"from first principles"', like the one we just finished, much more complicated. Instead of vanishing 
one-point moments, like (V£:(x) £:(x)) = 0, we would have to deal with correlation functions, like 
(Ve(x + r) e(x)), which generally do not vanish for r 7^ 0. As a result, we would not be able to 
reduce the dimensionality of the PDF through factorization formulae like equation (A14). In such 
cases it is much simpler to do as we did in section 3: use the Edgeworth expansion, with the values 
of cumulants (like S3) calculated from perturbation theory for the filtered field. 
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